disp('hqlfc');

Nc = 16;

vh = linspace(1e-4, 1e-1, 80);
vres = zeros(1, length(vh));

for ih = 1:length(vh)
    
    h0 = vh(ih);
    
    vk = linspace(0, 10, 100);
    mr = zeros(Nc, length(vk));
    
    for ik = 1:length(vk)
        
        k = vk(ik);
        mr(:,ik) = fn_solve_ft(k, Nc, h0);
        
    end
    
    % close(figure(1));
%     figure(1);
%     
%     for n = 1:Nc
%         vr = mr(n, :);
%         figure(1); hold on;
%         plot(vk, vr, '-o');
%     end
    
%     disp('lfll');
    
    mr2 = zeros(size(mr));
    
    for n = 1:Nc
        
        vr = mr(n, :);
        vplus = exp(-n/2*abs(vk));
        vr2 = vr.*vplus;
        mr2(n, :) = vr2*sinh(h0/2)^2/sinh(h0/2*(n+1))^2;
        
    end
    
    % close(figure(2));
    figure(2);
    
    for n = 1:Nc
        vr = mr2(n, :);
        figure(2); hold on;
        plot(vk, vr, '-o');
    end
    
    
    vr2 = sum(mr2);
    % figure(4); hold on;
    % plot(vk, tv, '-o');
    vk_nt = zeros(size(vk));
    vlen = vk(2:end) - vk(1:end-1);
    vk_nt(1:end-1) = 1/2*vlen;
    vk_nt(2:end) = vk_nt(2:end) + 1/2*vlen;
    
    % disp(sum(vk_nt.*vr2));
    vres(ih) = sum(vk_nt.*vr2);
    
end


% 
% 
% for ih = 1:length(vh)
%     
%     h0 = vh(ih);
%     
%     vk = linspace(0, 10, 800);
%     mr = zeros(Nc, length(vk));
%     
%     for ik = 1:length(vk)
%         
%         k = vk(ik);
%         mr(:,ik) = fn_solve_ft(k, Nc, h0);
%         
%     end
%     
%     
%     % close(figure(1));
%     % figure(1);
% %     for n = 1:Nc
% %         
% %         figure(1); hold on;
% %         plot(vk, mr(n,:), '-o');
% %         
% %     end
%     
%     mr2 = zeros(size(mr));
%     
%     for n = 1:Nc
%         
%         vr = mr(n, :);
%         vplus = exp(-n/2*abs(vk));
%         vr2 = vr.*vplus;
%         mr2(n, :) = vr2;
%         
%     end
%     
% %     for n = 1:Nc
% %         
% %         figure(2); hold on;
% %         plot(vk, mr2(n,:), '-o');
% %         
% %     end
%     
%     vr2 = sum(mr2);
%     
%     figure(3); hold on; 
%     plot(vk, vr2, '-o');
%     
%     vk_nt = zeros(size(vk));
%     vlen = vk(2:end) - vk(1:end-1);
%     vk_nt(1:end-1) = 1/2*vlen;
%     vk_nt(2:end) = vk_nt(2:end) + 1/2*vlen;
%     
%     disp(sum(vk_nt.*vr2));
%     
%     vres(ih) = sum(vk_nt.*vr2);
%     
% end
% 
% figure(999); hold on;
% plot(vh, vres, '-o');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function res = fn_solve_ft(k, Nc, h0)

Vn = fn_Vn(k, Nc);
An = fn_An(h0, Nc);
T = fn_T(k, Nc);

M = eye(Nc) + T*diag(An);

fprintf('cond(M) = %f \n', cond(M));

res = M\Vn;

end


function Vn = fn_Vn(k, Nc)

Vn = zeros(Nc, 1);

for n = 1:Nc
    
    Vn(n) = exp(-n/2*abs(k));
    
end

end


function An = fn_An(h0, Nc)

An = zeros(1, Nc);

for n = 1:Nc
    
    An(n) = sinh(h0/2)^2/sinh(h0/2*(n+1))^2;
    
end

end


function T = fn_T(k, Nc)

T = zeros(Nc, Nc);

if abs(k)<1e-10
    
    for n = 1:Nc
        for m = 1:Nc
            T(n,m) = 2*min(n,m);
        end
    end
    
else
    
    for n = 1:Nc
        for m = 1:Nc
            T(n,m) = (1+exp(-abs(k)))/(1-exp(-abs(k))) ...
                *(exp((-abs(n-m))/2*abs(k)) ...
                - exp((-abs(n+m))/2*abs(k)));
        end
    end
    
end

end